Set up libraries
Load data
c14dates <- read_csv("c14dates.csv")
Rows: 1452 Columns: 16── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (10): industry, site, unit, calib.curve, sample.id, material, delta 13, method&id, comments, citation
dbl (6): ID, Use, BP.cal.median, C14.mean, C14.SD, C14.COV
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rr.papers %>%
ggplot(aes(x=years)) +
geom_path(aes(y=JAS.risk/JAS.papers, color='Journal of Archaeological Science', lty='risk'), group=1, lwd=1.5) +
geom_path(aes(y=JAS.resilience/JAS.papers, color='Journal of Archaeological Science', lty='resilience'), group=1, lwd=1.5) +
geom_path(aes(y=Antiquity.risk/Antiquity.papers, color='Antiquity', lty='risk'), group=1, lwd=1.5) +
geom_path(aes(y=Antiquity.resilience/Antiquity.papers, color='Antiquity', lty='resilience'), group=1, lwd=1.5) +
scale_y_continuous(labels = scales::percent) +
scale_color_manual(values = c('Journal of Archaeological Science' = 'red', 'Antiquity' = 'blue')) +
scale_linetype_manual(values = c('risk' = 'solid', 'resilience' = 'dashed')) +
labs(title = "Papers Mentioning Risk or Resilience",
subtitle = "Change in Past 40 Years",
x="5-year intervals",
y="% of papers published",
color='journal:',
linetype='topic:') +
theme_bw(base_size = 30) +
theme(legend.position = 'bottom', legend.key.width = unit(50, "points"))
icecores %>%
ggplot() +
geom_point(aes(x=years.BP, y=d18O.NGRIP2.ppt), color='blue', size=.5, alpha=.2) +
geom_point(aes(x=years.BP, y=d18O.GISP2.ppt), color='pink2', size=.5, alpha=.2) +
geom_smooth(aes(x=years.BP, y=d18O.NGRIP2.ppt, color='ngrip2', lty='ngrip2'), method='loess', span=.01, se=FALSE, lwd = 1) +
geom_smooth(aes(x=years.BP, y=d18O.GISP2.ppt, color='gisp2', lty='gisp2'), method='loess', span=.01, se=FALSE, lwd = 1) +
scale_colour_manual(name=NULL, values =c('gisp2'='indianred3','ngrip2'='navy'), labels = c('GISP2','NGRIP2')) +
scale_linetype_manual(name=NULL, values =c('gisp2'='solid','ngrip2'='solid'), labels = c('GISP2','NGRIP2')) +
scale_x_reverse() +
#scale_x_continuous(breaks=c(19000,14000,10000), labels = c("19000","14000","10000"), trans = "reverse") +
geom_vline(xintercept = c(27000,22000,11000), lty='dashed', lwd=1) +
labs(x="years cal BP",
y="delta O18\ncolder << — >> warmer",
color="Cores",
title="Paleoclimate Proxies from Greenland Ice Cores") +
theme_minimal(base_size = 18) +
theme(legend.position="bottom",
plot.title = element_text(face = "bold", hjust = 0.5),
)
Hominin_demography %>%
mutate(fitness.aggregated = Fitness,
fitness.aggregated = recode(fitness.aggregated,
"birthrate=deathrate"="equal fitness",
"birthrate>deathrate"="equal fitness",
"MM higher mortality"="MM less fit",
"MM lower fertility"="MM less fit",
"MM lower mortality"="MM more fit",
"MN higher mortality"="MN less fit",
"MN lower fertility"="MN less fit",
"MN lower mortality"="MN more fit",
"NN higher mortality"="NN less fit",
"NN lower fertility"="NN less fit",
"NN lower mortality"="NN more fit"),
fitness.aggregated = factor(fitness.aggregated,
levels = c(
"equal fitness",
"MM less fit",
"MM more fit",
"MN less fit",
"MN more fit",
"NN less fit",
"NN more fit"))) %>%
ggplot(aes(x=phenotype, y=population)) +
geom_boxplot(fill="grey") +
facet_grid(fitness.aggregated~factor(home.range.radius), scales="free_y") +
labs(x="agent phenotype", y="population at end of simulation") +
theme_bw(base_size = 16)
Hominin_demography_trajectory %>%
ggplot(aes(y=population, x=step)) +
geom_line(aes(color=phenotype), size=2) +
scale_color_manual(values = c("red", "orange", "grey", "forestgreen", "blue")) +
labs(x="model step", y="population", color="agent phenotype") +
theme_bw(base_size = 16) +
theme(legend.position="bottom")
c14dates$BP.cal.median <- c14dates %>%
with(., calibrate(x=C14.mean, errors=C14.SD, calCurves = calib.curve, normalised=TRUE, calMatrix=FALSE)) %>%
medCal()
Dates file filtered to remove all date with COV > 0.1
par(mar=c(7,10,7,3))
plot(model.all, xlim = c(35000,8000), col.obs = 'red', lwd.obs = 3, drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("West Mediterranean Late Pleistocene Assemblages (N = ", nrow(subset(c14dates, calib.curve != "normal")), ")"), cex.main = 1.5)
c14dates <- c14dates %>%
mutate(industry.group =
case_when(
str_detect(tolower(industry), "aurig") ~ "Aurignacian",
str_detect(tolower(industry), "gdrav") ~ "Gravettian",
str_detect(tolower(industry), "solut") ~ "Solutrean",
str_detect(tolower(industry), "middle solutrean") ~ "Solutrean",
str_detect(tolower(industry), "salp") ~ "Solutrean",
str_detect(tolower(industry), "epig") ~ "Magdalenian",
str_detect(tolower(industry), "epi-gr") ~ "Magdalenian",
str_detect(tolower(industry), "eppig") ~ "Magdalenian",
str_detect(tolower(industry), "eppgr") ~ "Magdalenian",
str_detect(tolower(industry), "apig") ~ "Magdalenian",
str_detect(tolower(industry), "magd") ~ "Magdalenian",
str_detect(tolower(industry), "badeg") ~ "Magdalenian",
str_detect(tolower(industry), "gravet") ~ "Gravettian",
str_detect(tolower(industry), "gravt") ~ "Gravettian",
str_detect(tolower(industry), "azil") ~ "Epipaleolithic",
str_detect(tolower(industry), "romanell") ~ "Epipaleolithic",
str_detect(tolower(industry), "epipal") ~ "Epipaleolithic",
str_detect(tolower(industry), "mesol") ~ "Mesolithic",
str_detect(tolower(industry), "sauve") ~ "Mesolithic",
str_detect(tolower(industry), "tarden") ~ "Mesolithic",
str_detect(tolower(industry), "montcl") ~ "Mesolithic",
str_detect(tolower(industry), "montad") ~ "Mesolithic",
str_detect(tolower(industry), "castel") ~ "Mesolithic"
),
industry.group = factor(industry.group,
levels = c("Aurignacian",
"Gravettian",
"Solutrean",
"Magdalenian",
"Epipaleolithic",
"Mesolithic")),
.after=industry)
summary(perm.test)
'permTest()' function summary:
Number of sets: 6
Statistical Significance computed using 200 simulations.
--- Aurignacian ---
Number of radiocarbon dates:94
Number of bins:84
Global p-value: 0.00498
Significant positive local deviations at:
40000~29265 BP
28924~28146 BP
27436~26419 BP
Significant negative local deviations at:
18435~12504 BP
12048~11097 BP
10230~7646 BP
--- Gravettian ---
Number of radiocarbon dates:131
Number of bins:109
Global p-value: 0.00498
Significant positive local deviations at:
32846~25171 BP
24609~24371 BP
24318~24295 BP
24181~24163 BP
Significant negative local deviations at:
21508~21288 BP
20184~19623 BP
16356~13859 BP
13809~10550 BP
9793~7266 BP
--- Solutrean ---
Number of radiocarbon dates:120
Number of bins:88
Global p-value: 0.00498
Significant positive local deviations at:
26079~19630 BP
Significant negative local deviations at:
31656~31025 BP
18570~18569 BP
18560~18229 BP
18121~18117 BP
13728~12489 BP
12168~7419 BP
--- Magdalenian ---
Number of radiocarbon dates:516
Number of bins:350
Global p-value: 0.00498
Significant positive local deviations at:
20402~20034 BP
19782~12264 BP
Significant negative local deviations at:
36464~24702 BP
24446~23560 BP
23394~23225 BP
23221~22834 BP
11441~10975 BP
10860~7250 BP
--- Epipaleolithic ---
Number of radiocarbon dates:120
Number of bins:94
Global p-value: 0.00498
Significant positive local deviations at:
13305~10489 BP
10483~10341 BP
10146~9669 BP
Significant negative local deviations at:
32202~26954 BP
26507~25629 BP
25500~22554 BP
20524~19706 BP
18871~17393 BP
16130~14915 BP
--- Mesolithic ---
Number of radiocarbon dates:312
Number of bins:206
Global p-value: 0.00498
Significant positive local deviations at:
11324~7000 BP
Significant negative local deviations at:
34855~11918 BP
plot(perm.test,
focalm = "Aurignacian",
xlim = c(35000,8000),
col.obs = 'red',
lwd.obs = 3,
bbty = 'b',
ylim = c(0, 0.04),
drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("Aurignacian (N = ", nrow(subset(c14dates, industry.group == "Aurignacian")), ")"), cex.main = 1.5)
plot(perm.test,
focalm = "Gravettian",
xlim = c(35000,8000),
col.obs = 'red',
lwd.obs = 3,
bbty = 'b',
ylim = c(0, 0.04),
drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("Gravettian (N = ", nrow(subset(c14dates, industry.group == "Gravettian")), ")"), cex.main = 1.5)
plot(perm.test,
focalm = "Solutrean",
xlim = c(35000,8000),
col.obs = 'red',
lwd.obs = 3,
bbty = 'b',
ylim = c(0, 0.04),
drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("Solutrean (N = ", nrow(subset(c14dates, industry.group == "Solutrean")), ")"), cex.main = 1.5)
plot(perm.test,
focalm = "Magdalenian",
xlim = c(35000,8000),
col.obs = 'red',
lwd.obs = 3,
bbty = 'b',
ylim = c(0, 0.1),
drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("Magdalenian (N = ", nrow(subset(c14dates, industry.group == "Magdalenian")), ")"), cex.main = 1.5)
plot(perm.test,
focalm = "Epipaleolithic",
xlim = c(35000,8000),
col.obs = 'red',
lwd.obs = 3,
bbty = 'b',
ylim = c(0, 0.04),
drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("Epipaleolithic (N = ", nrow(subset(c14dates, industry.group == "Epipaleolithic")), ")"), cex.main = 1.5)
plot(perm.test,
focalm = "Mesolithic",
xlim = c(35000,8000),
col.obs = 'red',
lwd.obs = 3,
bbty = 'b',
ylim = c(0, 0.12),
drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)
axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)
mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("Mesolithic (N = ", nrow(subset(c14dates, industry.group == "Mesolithic")), ")"), cex.main = 1.5)